ELECTRO_ION

Overview

The ELECTRO_ION function fits specialized electrophysiology models to experimental data using non-linear least squares optimization. It supports four biophysically-inspired models commonly used in neuroscience and ion channel research: Boltzmann ion channel current-voltage relationships, synaptic current dynamics with triple rise and double decay components, action potential waveforms, and the Goldman-Hodgkin-Katz (GHK) membrane potential equation.

The function leverages scipy.optimize.curve_fit from SciPy to estimate optimal parameters for each model. For the Boltzmann ion channel current-voltage model, the relationship between membrane voltage and current is described by:

I(V) = \frac{(V - V_{rev}) \cdot g_{max}}{1 + e^{(V - V_{half})/dx}}

where V_{half} is the half-activation voltage, dx is the slope factor, g_{max} is maximum conductance, and V_{rev} is the reversal potential. This sigmoid-shaped function models the voltage-dependent activation of ion channels that underlies neuronal excitability.

The Goldman-Hodgkin-Katz membrane potential model calculates the equilibrium potential across a cell membrane based on ionic concentrations and permeabilities:

E_m = \frac{RT}{F} \ln\left(\frac{[K^+]_o + b[Na^+]_o}{[K^+]_i + b[Na^+]_i}\right)

where b is the relative sodium-to-potassium permeability ratio, and R, T, and F are the gas constant, temperature, and Faraday constant respectively. This foundational equation, developed by David Goldman and later refined by Hodgkin and Katz, describes how the resting membrane potential emerges from differential ion permeabilities.

The synaptic current and action potential models capture complex temporal dynamics using exponential rise and decay functions with piecewise definitions before and after peak amplitude times. These models are essential for quantifying parameters from electrophysiological recordings such as patch clamp experiments.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=ELECTRO_ION(xdata, ydata, electro_ion_model)

xdata (list[list], required): The xdata value
ydata (list[list], required): The ydata value
electro_ion_model (str, required): The electro_ion_model value

Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error string.

Examples

Example 1: Demo case 1

Inputs:

electro_ion_model	xdata	ydata
boltzmann_ion_channel_current_voltage	0.01	-8.667462978597941
	2.0075	-4.517640100090819
	4.005	-1.1317796397731288
	6.0024999999999995	1.2016629480589822
	8	2.073117520844434

Excel formula:

=ELECTRO_ION("boltzmann_ion_channel_current_voltage", {0.01;2.0075;4.005;6.0024999999999995;8}, {-8.667462978597941;-4.517640100090819;-1.1317796397731288;1.2016629480589822;2.073117520844434})

Expected output:

vhalf	dx	gmax	vrev
5.895	2.848	2.027	4.822
0.2831	0.2666	0.0651	0.03723

Example 2: Demo case 2

Inputs:

electro_ion_model	xdata	ydata
synaptic_current_triple_rise_double_decay	0.01	3.135937877328851
	2.0075	5.301859757942096
	4.005	3.569379720613383
	6.0024999999999995	2.320611288976838
	8	1.6205926862209843

Excel formula:

=ELECTRO_ION("synaptic_current_triple_rise_double_decay", {0.01;2.0075;4.005;6.0024999999999995;8}, {3.135937877328851;5.301859757942096;3.569379720613383;2.320611288976838;1.6205926862209843})

Expected output:

y0	xc	Ag1	tg1	Ag2	tg2	Ag3	tg3	Ad1	td1	Ad2	td2
1.263	4.0259	1.1116	0.0679	0.384	0.0825	0	28.0457	1.2805	1.5022	1.8922	2.0285

Example 3: Demo case 3

Inputs:

electro_ion_model	xdata	ydata
action_potential_single_rise_decay	0.01	1.9684261274918025
	2.0075	2.6519549335611137
	4.005	1.785374148442621
	6.0024999999999995	1.159640232264233
	8	0.8106394055832311

Excel formula:

=ELECTRO_ION("action_potential_single_rise_decay", {0.01;2.0075;4.005;6.0024999999999995;8}, {1.9684261274918025;2.6519549335611137;1.785374148442621;1.159640232264233;0.8106394055832311})

Expected output:

y0	xc	Ag	tg	Ad	td
0.3705	2.37	0.7228	0.1792	2.281	3.421

Example 4: Demo case 4

Inputs:

electro_ion_model	xdata	ydata
goldman_hodgkin_katz_membrane_potential	0.1	0
	1.3250000000000002	0
	2.5500000000000003	0
	3.7750000000000004	0
	5	0

Excel formula:

=ELECTRO_ION("goldman_hodgkin_katz_membrane_potential", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0;0;0;0;0})

Expected output:

b	Ko	Nao	Ki	Nai	T
1.052	4.999	139.8	140.8	10.81	292.7

Python Code

import numpy as np
from scipy.optimize import curve_fit as scipy_curve_fit
import math

def electro_ion(xdata, ydata, electro_ion_model):
    """
    Fits electro_ion models to data using scipy.optimize.curve_fit. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html for details.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        xdata (list[list]): The xdata value
        ydata (list[list]): The ydata value
        electro_ion_model (str): The electro_ion_model value Valid options: Boltzmann Ion Channel Current Voltage, Synaptic Current Triple Rise Double Decay, Action Potential Single Rise Decay, Goldman Hodgkin Katz Membrane Potential.

    Returns:
        list[list]: 2D list [param_names, fitted_values, std_errors], or error string.
    """
    def _validate_data(xdata, ydata):
        """Validate and convert both xdata and ydata to numpy arrays."""
        for name, arg in [("xdata", xdata), ("ydata", ydata)]:
            if not isinstance(arg, list) or len(arg) < 2:
                raise ValueError(f"{name}: must be a 2D list with at least two rows")
            vals = []
            for i, row in enumerate(arg):
                if not isinstance(row, list) or len(row) == 0:
                    raise ValueError(f"{name} row {i}: must be a non-empty list")
                try:
                    vals.append(float(row[0]))
                except Exception:
                    raise ValueError(f"{name} row {i}: non-numeric value")
            if name == "xdata":
                x_arr = np.asarray(vals, dtype=np.float64)
            else:
                y_arr = np.asarray(vals, dtype=np.float64)

        if x_arr.shape[0] != y_arr.shape[0]:
            raise ValueError("xdata and ydata must have the same number of rows")
        return x_arr, y_arr

    # Model definitions dictionary
    models = {
        'boltzmann_ion_channel_current_voltage': {
            'params': ['vhalf', 'dx', 'gmax', 'vrev'],
            'model': lambda x, vhalf, dx, gmax, vrev: ((x - vrev) * gmax) / (1.0 + np.exp((x - vhalf) / dx)),
            'guess': lambda xa, ya: (float(np.median(xa)), 1.0, float(np.ptp(ya) / (np.ptp(xa) if np.ptp(xa) else 1.0)), float(np.min(xa))),
            'bounds': ([-np.inf, -np.inf, 0.0, -np.inf], np.inf),
        },
        'synaptic_current_triple_rise_double_decay': {
            'params': ['y0', 'xc', 'Ag1', 'tg1', 'Ag2', 'tg2', 'Ag3', 'tg3', 'Ad1', 'td1', 'Ad2', 'td2'],
            'model': lambda x, y0, xc, Ag1, tg1, Ag2, tg2, Ag3, tg3, Ad1, td1, Ad2, td2: np.where(x <= xc, y0 + Ad1 + Ad2 + Ag1 * (np.exp(-xc / tg1) - np.exp(-x / tg1)) + Ag2 * (np.exp(-xc / tg2) - np.exp(-x / tg2)) + Ag3 * (np.exp(-xc / tg3) - np.exp(-x / tg3)), y0 + Ad1 * np.exp(-(x - xc) / td1) + Ad2 * np.exp(-(x - xc) / td2)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), 1.0, 1.0, 0.5, 2.0, 0.5, 3.0, 0.5, 1.0, 0.5, 2.0),
            'bounds': ([-np.inf, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], np.inf),
        },
        'action_potential_single_rise_decay': {
            'params': ['y0', 'xc', 'Ag', 'tg', 'Ad', 'td'],
            'model': lambda x, y0, xc, Ag, tg, Ad, td: np.where(x <= xc, y0 + Ad + Ag * (np.exp(-xc / tg) - np.exp(-x / tg)), y0 + Ad * np.exp(-(x - xc) / td)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), 1.0, 1.0, 1.0, 1.0),
            'bounds': ([-np.inf, 0.0, 0.0, 0.0, 0.0, 0.0], np.inf),
        },
        'goldman_hodgkin_katz_membrane_potential': {
            'params': ['b', 'Ko', 'Nao', 'Ki', 'Nai', 'T'],
            'model': lambda x, b, Ko, Nao, Ki, Nai, T: np.full_like(x, (8.314462618 * T / 96485.33212) * np.log((Ko + b * Nao) / (Ki + b * Nai))),
            'guess': lambda xa, ya: (0.5, 5.0, 140.0, 140.0, 10.0, 293.0),
            'bounds': (0.0, np.inf),
        }
    }

    # Validate model parameter
    if electro_ion_model not in models:
        return f"Invalid model: {str(electro_ion_model)}. Valid models are: {', '.join(models.keys())}"

    model_info = models[electro_ion_model]

    # Validate and convert input data
    try:
        x_arr, y_arr = _validate_data(xdata, ydata)
    except ValueError as e:
        return f"Invalid input: {e}"

    # Perform curve fitting
    try:
        p0 = model_info['guess'](x_arr, y_arr)
        bounds = model_info.get('bounds', (-np.inf, np.inf))
        if bounds == (-np.inf, np.inf):
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, maxfev=10000)
        else:
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, bounds=bounds, maxfev=10000)

        fitted_vals = [float(v) for v in popt]
        for v in fitted_vals:
            if math.isnan(v) or math.isinf(v):
                return "Fitting produced invalid numeric values (NaN or inf)."
    except ValueError as e:
        return f"Initial guess error: {e}"
    except Exception as e:
        return f"curve_fit error: {e}"

    # Calculate standard errors
    std_errors = None
    try:
        if pcov is not None and np.isfinite(pcov).all():
            std_errors = [float(v) for v in np.sqrt(np.diag(pcov))]
    except Exception:
        pass

    return [model_info['params'], fitted_vals, std_errors] if std_errors else [model_info['params'], fitted_vals]

Overview

Excel Usage

Examples

Python Code

Online Calculator